# ===============================================================================================================================
# ===============================================================================================================================
# Clarity Analysis:
# This code is for replicating: Figures 4 and 5, Tables 3, 4, and 5,  Appendix Tables A7, A8, A9, A10, A11, and C1
# 
# This code runs all of the Clarity analysis using plasma opening cohorts.
# The code is split into several groupings (largely based on the dataset being used):
#   1. CIFP cohort inquiry full panel
#   2. CTFP cohort transactions full panel
#   3. CTFE cohort transactions event anlysis (collapsed to simple pre-post analysis)
#   4. Create tables 
#   5. Create figures
# ===============================================================================================================================
# ===============================================================================================================================


# ===============================================================================================================================
# ===============================================================================================================================
# Establish directories
# ===============================================================================================================================
# ===============================================================================================================================
#Clear out the namespace
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects. (simpler is rm(list=ls()))
gc() #free up memory and report the memory usage.
.pardefault <- par()

#Turn off scientific notation (http://datacornering.com/remove-scientific-notation-in-r/)
# options(scipen = 100)
options(scipen = 999)

# Set the path to the root folder ending in RFSDelivery 
cd <- "D:/..."
# Initialize all other paths and import all functions and packages from the AllSettings file.
source(paste0(cd,'/Code/R/AllSettings.R',collapse=''))

#Turn graphics off to allow saving the plots
graphics.off()
#List of coefficients which are estimated and included in the plots 
dyncoef_drop = "(-16|-15|-14|-13|-12|-11|-10|-9|13|14|15)"

# version control
cfolder <- ''
scriptv <- 'v8'

# ===============================================================================================================================
# ===============================================================================================================================
# Import the Clarity data set
# ===============================================================================================================================
# ===============================================================================================================================
# Import the variable dictionary for labeling heterogeneity variables
cla_desc <- read_excel(cdd['p_d_cl'] %c% '/ClarityVariables.xlsx', sheet='ValLabels')
cla_desc <- cla_desc %>% mutate(value = as.numeric(value))
labels_het <- cla_desc %>% distinct(Alias, description, .keep_all = TRUE)
labels_het <- as.list(setNames(labels_het$description,labels_het$variable))

################Import the banking desert data###################
# The "banking desert" data (commdata_revised.dta) was provided by Friedline, Despard, and West (2017). Bank branch locations were collected through the FDIC's summary
# of deposits. Credit union branch locations were collected through the NCUA call reports. We use zip code-level aggregates of traditional banks and
# credit unions according to regulatory filings as of 2014.
bd <- read_dta(cdd['p_d'] %c% '/Banking Deserts/commdata_revised.dta') 
bd <- bd %>% rename(zip = zipcode) %>% mutate(desert = replace_na(desert, 0),
                                              altfin = (replace_na(TotalAFS2, 0) - replace_na(TaxPrep2 , 0) - replace_na(MoneyOrder2, 0)) / TOTPOP_CY ,
                                              num_altfin = (replace_na(TotalAFS2, 0) - replace_na(TaxPrep2 , 0) - replace_na(MoneyOrder2, 0))  ,
                                              unemp = UNEMP_CY / TOTPOP_CY,
                                              nonwhite = 1 - (WHITE_CY/ TOTPOP_CY),
                                              black = (BLACK_CY/ TOTPOP_CY), 
                                              banked = (replace_na(NCUA2014, 0) + replace_na(FDIC2014, 0)) / TOTPOP_CY,
                                              num_banks = (replace_na(NCUA2014, 0) + replace_na(FDIC2014, 0))
                                              )
bd <- bd %>% select(zip,desert,altfin,num_altfin,TotalAFS2,TaxPrep2,TOTPOP_CY,MoneyOrder2,unemp,UNEMP_CY,UNEMPRT_CY,nonwhite,WHITE_CY,black,BLACK_CY,banked,NCUA2014,FDIC2014,MEDHINC_CY,num_banks) 


# Variable labels
labels_credit_utilization = list(
                                  st_gc_f='Loans Good',st_del_f='Loans Del',st_col_f='Loans Col',st_co_f='Loans CO'
                                )
labels_inq = list(
                    inq_payday_i_any='Payday Int gt0',inq_inst_i_any='Installment Int gt0',
                    inq_payday_s_any='Payday Phys gt0',inq_inst_s_any='Installment Phys gt0',
                    inq_payday_any='Payday gt0',inq_inst_any='Installment gt0',
                    inq_payday_gt1='Payday gt1',inq_inst_gt1='Installment gt1',
                  )
labels_trans = list(
                    ntrans_da_payday_any='Payday gt0',
                  )
labels_transe = list(
                      bal_highprin_w_mu='Payday Bal',bal_highprin_w_mx='Payday Bal',
                      paidlate_mu='Late Pmt Rate',chargeoff_mx='Charge-Off Rate'
                    )
labels_coh = list(treated='Treated',post='Post',etime='Event Time',etimey='Event Time',distance_t2cf='Dist of Target PC to CF PC',
                  date='Date',ctype='Counterfactual type',cohort='Cohort')
labels_indiv = list(id='Indiv',zip='Zip Code',state='State',housing='Housing',pay_freq='Pay Freq',
                    dor='Date of Residence',dob='Date of Birth',
                    age='Age',reside='Time at Residence',income='Net Monthly Income',moved='Moved',
                    age_2g='Age lte35',income_2g='Income lte 2K',age_3g='Age 3Q',income_3g='Income 3Q')
labels_indep = list(distance='Dist',ic='IC')
labels_treat = list(d50_ic50='(D<5)(IC>5)',na='All')
labels_other = list(pop2sqkm_wq10='Density')

labels_all <- c(labels_credit_utilization, labels_inq, labels_inqe, labels_trans, labels_transe, 
                labels_coh, labels_indiv, labels_indep, labels_other, labels_treat)


#Set the labels used for all tables
my.style.df = style.df(fixef.title = "", fixef.suffix = " FE", stats.title = "", slopes.title = "")
my.style.tex = style.tex(main="base", fixef.suffix = " FE", slopes.title = "",
                         tablefoot = FALSE,line.bottom="\\bottomrule",line.top="\\toprule",
                         depvar.title = "",var.title="\\midrule", fixef.title = '\\midrule', stats.title = '\\midrule')
setFixest_etable(digits = "r4", digits.stats = "r4", style.df = my.style.df, style.tex = my.style.tex, fitstat = c('n','r2'))
setFixest_dict(unlist(labels_all))



# ===============================================================================================================================
# ===============================================================================================================================
# SECTION 1: INQUIRIES ANALYSIS
# ===============================================================================================================================
# ===============================================================================================================================

###########################Begin Inquiry dataset build###########################
cifp <- read_parquet(cdd['p_d_cl_pa'] %c% cfolder %c% '/cifp_q.parquet')
t <- toc()

# Filter the data
cifp <- cifp %>% group_by(cohort) %>% mutate(anytreated = max(treated)) %>% ungroup
cifp <- cifp %>% filter(anytreated == 1)
# #Restrict observations based on the available data
cifp <- cifp %>% filter((open >= as.Date('2014-07-01'))&(open < as.Date('2020-01-01')))
cifp <- cifp %>% filter((distance<5000)&(ic>5000))


# I choose to fill all na with zeros. The only column with NA is C13  which has a couple years missing. But since I study this alongside C5 it should be fine for DiD analysis.
cifp <- cifp %>% mutate_at(vars(names(cifp)[grep('inq',names(cifp))]),~replace(., is.na(.), 0))  %>%
           mutate(
                      inq_payday = inq_C1 + inq_C2,
                      inq_inst = inq_C3 + inq_C6,
                      inq_payday_any = 1*(inq_payday>0),
                      inq_inst_any = 1*(inq_inst>0),
                      inq_payday_i_any = 1*(inq_C1>0),
                      inq_payday_s_any = 1*(inq_C2>0),
                      inq_inst_i_any = 1*(inq_C3>0),
                      inq_inst_s_any = 1*(inq_C6>0),
                      d50_ic50_19=1*((distance<5000)&(ic>5000)&(date<as.Date('2020-01-01'))&(open<as.Date('2020-01-01'))),
                      etimey = floor(etime /4),
                      post = 1*(etime>=0),
                      na=1
                    )

cifp <- cifp %>% mutate(
                        age_2gr = 1-age_2g,
                        ct_past = 1*((ctype==5)|(ctype==2)),
                        ct_future = 1*((ctype==5)|(ctype==1)),
                        ct_fp2 = 1*((ctype==5)|(ctype==1)|(difftime(open ,open_cf , units = c("days"))>8*365))
                        )


##################Merge in banking desert data##########################
cifp <- left_join(cifp,bd %>% select(zip,desert,altfin,unemp,UNEMPRT_CY,nonwhite,black,banked,MEDHINC_CY,num_altfin,num_banks),by='zip')

# Create some factor variables.
cifp <- cifp %>% mutate(post = 1*(etime>=0), na=1)
cifp <- cifp %>% group_by(cohort,id,post) %>% mutate(
                                                    inq_payday_npre = sum(inq_payday),
                                                     inq_inst_npre = sum(inq_inst)
                                                     ) %>% ungroup()
cifp <- cifp %>% arrange(cohort,id,etime) %>% mutate(
                          inq_payday_npre = case_when(post==0 ~ inq_payday_npre, post==1 ~ NA_real_),
                          inq_inst_npre = case_when(post==0 ~ inq_inst_npre, post==1 ~ NA_real_),
                          ) %>% group_by(cohort,id) %>%  fill(inq_payday_npre,inq_inst_npre) %>% ungroup()


###################heterogeneity split variables########################
# cut split variables into binary variables
cifp <- cifp %>% mutate(
                          inq_payday_npre_any = if_else(is.na(inq_payday_npre>0),0,1*(inq_payday_npre>0)),
                          inq_inst_npre_any = if_else(is.na(inq_inst_npre>0),0,1*(inq_inst_npre>0)),
                          inq_npre_any = 1*(inq_payday_npre_any+inq_inst_npre_any>0),
                          inq_npre_gt1 = 1*(inq_payday_npre+inq_inst_npre>1),
                          inq_pd_npre_any = 1*(inq_payday_npre_any>0),
                          inq_instal_npre_any = 1*(inq_inst_npre_any>0),

                          unemprt_2g = 1*(UNEMPRT_CY > median(cifp[['UNEMPRT_CY']],na.rm=TRUE)),
                          unemprt_bot4g = 1*(UNEMPRT_CY < 5.5),
                          unemprt_top4g = 1*(UNEMPRT_CY > 10.1),
                          medinc_2g = 1*(MEDHINC_CY > median(cifp[['MEDHINC_CY']],na.rm=TRUE)),
                          medinc_bot4g = 1*(MEDHINC_CY < 33891),
                          medinc_top4g = 1*(MEDHINC_CY > 53703),
                          nonwhite_2g = 1*(nonwhite > median(cifp[['nonwhite']],na.rm=TRUE)),
                          nonwhite_bot4g = 1*(nonwhite < 0.27568),
                          nonwhite_top4g = 1*(nonwhite > 0.65049),
                          black_2g = 1*(black > median(cifp[['black']],na.rm=TRUE)),
                          black_bot4g = 1*(black < 0.05939977),
                          black_top4g = 1*(black > 0.42227012),
                          altfin_2g = 1*(altfin > median(cifp[['altfin']],na.rm=TRUE)),
                          altfin_bot4g = 1*(altfin < 0.0001097775),
                          altfin_top4g = 1*(altfin > 0.0003016981),
                          banked_2g_lo = 1*(banked < median(cifp[['banked']],na.rm=TRUE)),
                          banked_bot4g = 1*(banked < 0.0001383158),
                          banked_top4g = 1*(banked > 0.0003936620),
                          
                          num_altfin_2g = 1*(num_altfin > median(cifp[['num_altfin']],na.rm=TRUE)),
                          num_altfin_bot4g = 1*(num_altfin < 4),
                          num_altfin_top4g = 1*(num_altfin > 13),
                          num_banked_2g = 1*(num_banks > median(cifp[['num_banks']],na.rm=TRUE)),
                          num_banked_bot4g = 1*(num_banks < 5),
                          num_banked_top4g = 1*(num_banks > 15),
                          
                          density_urban = 1*(pop2sqkm_wq3==2),
                          density_topurban = 1*(pop2sqkm_wq10>7),
                          density_gt50urban = 1*(pop2sqkm_wq10>4)
                        )

# Apply variable labels
labels_sub <- labels_all[names(labels_all) %in% names(cifp)]
var_label(cifp) <- labels_sub


###################set up regressions we want to loop through########################## 
spec_task = c('est','tbl_sum')
#Spec_y is outcomes.	Note that there is no need to specify a set of outcome variables spec_y. fixest will estimate a joint regression for all of the desired outcomes at once. Even if there are several groups of outcomes specified this is purely for legibility. My code will concatenate all non-commented lists of outcomes into a single list for joint estimation.
spec_y = list( 
              # MAIN:
              'inq_simpa'=c('inq_payday_any','inq_inst_any'),
              # storefront v. internet TEST:
               'inq_simpg'=c('inq_payday_s_any','inq_payday_i_any','inq_inst_s_any','inq_inst_i_any')
            )
#Spec_i is independent variables
spec_i = list(
                'Tty' = 'i(etimey, treated, ref = -1)'
                # 'Tp' = 'treated:post'
              )
#Spec_fe is fixed effects
spec_fe = list(
                  'ci_ctpd10_sd'='cohort^id + cohort^etime^pop2sqkm_wq10 + state^date'
                )
#Dictionary to handle the splits by defining the variable and value labels and the base category (for marginal effects)
# Spec_split is the heterogeneity variables we are cutting on. Must edit according to the cuts you want to do
spec_split = list(
                   'na'='na',
                  # Figure 4 and Appendix Table A7
                   'age_2g_b1'=list('var'='age_2g','var_label'='Age','base'=1,'val_label'=c('Age lte35'=0,'Age gt35'=1)),
                  # Table 3
                    'income_3g_low'=list('var'='income_3g','var_label'='Income','base'=2,'val_label'=c('Inc lt2K'=0,'Inc 2-3K'=1,'Inc gt3K'=2)),
                  # Table A11
                    'density_3g_low'=list('var'='pop2sqkm_wq3','var_label'='Income','base'=2,'val_label'=c('Rural'=0,'Suburban'=1,'Urban'=2)),
                  # Table A10 
                    'num_inq_npre_gt1'=list('var'='inq_npre_gt1','var_label'='Inquiries','base'=0,'val_label'=c('Inq gt1'=0,'Inq gt1'=1)),
                    'banked_2g_low'=list('var'='banked_2g_lo','var_label'='Banks','base'=0,'val_label'=c('Banks High'=0,'Banks Low'=1)),
                    'unemprt_4g_top75'=list('var'='unemprt_top4g','var_label'='Unemp Rate','base'=0,'val_label'=c('Unemp Low'=0,'Unemp top75'=1)),
                  # Table C1
                    'yinq_npre_any'=list('var'='inq_npre_any','var_label'='Inquiries','base'=0,'val_label'=c('Inq eq0'=0,'Inq gt0'=1))
                )
#    e.g., the split variable will be banked_2g_lo where the base or omitted category is set to be when the variable 
#    banked_2g_lo takes the value 0 (high). 


# Spec_split_m is the method of testing for splits 
#   "s" is a regression for each subsample (a split sample)
#   "ig" and "is" test for incremental effect of each group compared to a baseline group value specified in the "spec_split" definition. 
#        "ig" creates a temporary generic variable name that can be swapped out later in tables 
#        "is" uses the exact variable name passed for tests making tables harder to generate later (so I recommend only using "ig").
spec_split_m = c('na','ig')  #  's','ig' 'na','s','is','ig'
# no need to bother with the other terms (interaction_types, spec_treat_rest, spec_cluster).
interaction_types = list('is'='specific','ig'='generic')
spec_treat_rest <- c(
                      # 'age_2gr,
                      # 'ct_past',
                      # 'ct_future',
                      # 'ct_fp2',
                      'd50_ic50_19'            #Full sample
                    )
spec_cluster <- list(
                      'coh'='cohort'
                    )
# spec_trl <- 'l'  #s for small: restricts to observations that are matched on treated and control fixed effect groups (versus l that includes everything)
spec_lean <- TRUE  #TRUE means small files or FALSE for large files.  (I HAVE VERIFIED THAT ETABLE USES THE CLUSTERING SPECIFIED AT ESTIMATION EVEN WITH LEAN)
spec_suf <- '_ds(cifp)'

gc()
mem_used()


##############run regressions specified above and save output################
for (var_tr in spec_treat_rest) {
  cifpt <- cifp %>% filter(get(var_tr)==1)
  for (var_s in names(spec_split)){
    for (var_spm in spec_split_m) {
      for (var_i in names(spec_i)){
        for (var_c in names(spec_cluster)){
          for (var_fe in names(spec_fe)){
            if ((var_s == 'na')&(var_spm!='na')){ next }
            if ((var_s != 'na')&(var_spm=='na')){ next }
            # Create folder for model, tables, figures
            specl <- ifelse(spec_lean, 'L', 'F')
            spec_ifes <- 'i(' %c% var_i %c% ')_fe('%c% var_fe %c% ')_c('%c% var_c %c%')_s('%c% var_s %c%')_sm(' %c% var_spm %c% ')' %c% '_trl(' %c% var_tr %c% ')_l(' %c% specl %c% ')_' %c% scriptv %c% spec_suf
            print(spec_ifes)
            mod <- 'DDiD_' %c% spec_ifes %c% '.Rqs'
            models <- list.files(path = cdd['p_rm'] %c% '/Clarity' %c% cfolder,pattern='\\.Rqs')
            filem <- cdd['p_rm_cl']  %c% cfolder %c% '/' %c% mod
            folder_rt <- cdd['p_rt'] %c% '/Clarity/' %c% cfolder %c% spec_ifes
            folder_rf <- cdd['p_rf'] %c% '/Clarity/' %c% cfolder %c% spec_ifes
            folder_rta <- cdd['p_rt'] %c% '/Clarity' %c% cfolder %c% '/All'
            folder_rfa <- cdd['p_rf'] %c% '/Clarity' %c% cfolder %c% '/All'
            
            var_y <- unlist(spec_y)
            var_y_missing <- var_y[!(var_y %in% names(cifp))]
            if (length(var_y_missing)>0){ 
              print(var_y_missing) 
              next
            }
            var_y <- var_y[var_y %in% names(cifp)]
            var_y <- 'c(' %c% paste0(var_y,collapse=', ') %c% ') '
            if ('est' %in% spec_task){
              res <- my_feols(
                df=cifpt,
                labels=labels_sub,
                form=list('y'=var_y,
                          'i'=spec_i[[var_i]],
                          'fe'=spec_fe[[var_fe]],
                          'cluster'=spec_cluster[[var_c]],
                          's'=spec_split[[var_s]],
                          'spm'=var_spm),
                filem=filem,
                lean=spec_lean
              )
              if (res$rc=='error'){
                print('estimation error')
                next
              }
              est <- res$est
            }
            if ('tbl_sum' %in% spec_task){
              tbl_all <- my_etable_df(est)
              openxlsx::write.xlsx(as_tibble(tbl_all), folder_rta %c% '/DDiD_' %c% spec_ifes %c% '.xlsx', row.names=FALSE, overwrite=TRUE)
            }
            gc()
            mem_used()
          }
        }
      }
    }
  }
}

# include these lines when you run the regression code above so that it will delete the massive datasets when the regressions are finished 
# and free up RAM for other users (to be considerate). That way you don't have to babysit the code.
rm(cifpt)
rm(cifp)
gc()
mem_used()

###########output naming convention##############
#	DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(density_3g_low)_sm(s)_trl(d50_ic50_19)_l(L)_v9_ds(cifp)_rnd4
#	DDiD is dynamic Diff-in-Diff
#	i(Tty) is independent variable is Tty which is defined in the code: 'Tty' = 'i(etimey, treated, ref = -1)'
#	fe(ci_ctpd10_sd) is the fixed effect specification defined:  'ci_ctpd10_sd'='cohort^id + cohort^etime^pop2sqkm_wq10 + state^date'
#	c(coh) is the clustering chosen: 'coh'='cohort'
#	s(d50_ic50_19) is the split variable
#	sm(ig) is the split method (i.e., split subsample or a test for a significant incremental treatment effect on a subpopulation).
#	Trl(age_2gr) is based on spec_treat_rest: 'd50_ic50_19' #Full sample.
#      It is restricting the sample to a subset of individuals (age_2gr is defined above and takes a value of 1 only for individuals younger than 35). So this would be testing for heterogeneity in the group with material treatment. I forget the exact shorthand I had in mind with trl but probably "treatment restriction list" or something thereabouts.
#	l(L) is an indicator for whether the regression saved a lean dataset. This is an option within fixest. By default it returns a model object with a copy of the dataset sufficient to change the clustering after estimation (effectively a full copy of the dataset that may be GB large). This is slow and takes a lot of space. So instead we always want to use the "lean" option that doesn't save a dataset with the model object. This also means you can't change clustering after estimation which is a cost worth paying.
#	Finally v8 is the version of the code and ds(cifp) is the dataset that you are studying.




# ===============================================================================================================================
# ===============================================================================================================================
# SECTION 2: TRANSACTIONS ANALYSIS
# ===============================================================================================================================
# ===============================================================================================================================

###########################Begin Transaction dataset build###########################
ctfp <- read_parquet(cdd['p_d_cl_pa'] %c% cfolder %c% '/ctfp_q.parquet')

# Filter the data
ctfp <- ctfp %>% group_by(cohort) %>% mutate(anytreated = max(treated)) %>% ungroup
ctfp <- ctfp %>% filter(anytreated == 1)
# #Restrict observations based on the available data
ctfp <- ctfp %>% filter((open >= as.Date('2014-07-01'))&(open < as.Date('2020-01-01')))
ctfp <- ctfp %>% filter((distance<5000)&(ic>5000))


# Create some factor variables.
ctfp <- ctfp %>% mutate(
                          ntrans_da_payday        = ntrans_daC1 + ntrans_daC2,
                          ntrans_da_inst          = ntrans_daC3 + ntrans_daC6,
                          ntrans_da_unsec         = ntrans_da01,
                          ntrans_da_payday_any    = 1*(ntrans_da_payday > 0),
                          ntrans_da_inst_any      = 1*(ntrans_da_inst > 0),
                          ntrans_da_unsec_any     = 1*(ntrans_da_unsec > 0),
                          ntrans_da_payday_i_any    = 1*(ntrans_daC1 > 0),
                          ntrans_da_payday_s_any    = 1*(ntrans_daC2 > 0),
                          ntrans_da_inst_i_any      = 1*(ntrans_daC3 > 0),
                          ntrans_da_inst_s_any      = 1*(ntrans_daC6 > 0),
                          d50_ic50_19=1*((distance<5000)&(ic>5000)&(date<as.Date('2020-01-01'))&(open<as.Date('2020-01-01'))),
                          etimey = floor(etime /4),
                          post = 1*(etime>=0),
                          na=1
                        )
# Create this variable to condition on individuals that have taken out a payday loan at least once.
ctfp <- ctfp %>% group_by(cohort, id) %>% mutate(ntrans_da_payday_anyc = 1*(sum(ntrans_da_payday)>0)) %>% ungroup()

# Must build a pre vs. post collapsed dataset for the regressions in Table 5
ctfp <- ctfp %>% group_by(cohort,id,post) %>% mutate(
                              payday_pre = 1*((sum(ntrans_da_payday)>0)&(post==0)),
                              payday_post = 1*((sum(ntrans_da_payday)>0)&(post==1))
                              ) %>% ungroup() %>% group_by(cohort,id) %>% mutate(
                              payday_pre = 1*(max(payday_pre)==1),
                              payday_post = 1*(max(payday_post)==1),
                              payday_prepost = 1*(payday_pre+payday_post==2)
                              ) %>% ungroup()

###################set up regressions we want to loop through########################## 

spec_task = c('est','tbl_sum')
spec_y = list(
                #Figure 5 and Appendix Table A8
                'ntrans_da_any'=c('ntrans_da_payday_any')
              )
spec_i = list(
                'Tty' = 'i(etimey, treated, ref = -1)'
                # 'Tp' = 'treated:post'
              )
spec_fe = list(
                'ci_ctpd10_sd'='cohort^id + cohort^etime^pop2sqkm_wq10 + state^date'
                )
#Dictionary to handle the splits by defining the variable and value labels and the base category (for marginal effects)
spec_split = list(
                    'age_2g_b1'=list('var'='age_2g','var_label'='Age','base'=1,'val_label'=c('Age lte35'=0,'Age gt35'=1)),
                    'na'='na'
                  )
spec_split_m = c('s','ig')  #'na','s','is','ig'
interaction_types = list('is'='specific','ig'='generic')
spec_treat_rest <- c(
                      'd50_ic50_19'            #Full sample
                    )
spec_cluster <- list(
                      'coh'='cohort'
                    )
# spec_trl <- 'l'  #s for small: restricts to observations that are matched on treated and control fixed effect groups (versus l that includes everything)
spec_lean <- TRUE  #TRUE means small files or FALSE for large files.  (I HAVE VERIFIED THAT ETABLE USES THE CLUSTERING SPECIFIED AT ESTIMATION EVEN WITH LEAN)
spec_suf <- '_ds(ctfp)'


##############run regressions specified above and save output################
for (var_tr in spec_treat_rest) {
  ctfpt <- ctfp %>% filter(get(var_tr)==1)
  for (var_s in names(spec_split)){
    for (var_spm in spec_split_m) {
      for (var_i in names(spec_i)){
        for (var_c in names(spec_cluster)){
          for (var_fe in names(spec_fe)){
            if ((var_s == 'na')&(var_spm!='na')){ next }
            if ((var_s != 'na')&(var_spm=='na')){ next }
            # Create folder for model, tables, figures
            specl <- ifelse(spec_lean, 'L', 'F')
            spec_ifes <- 'i(' %c% var_i %c% ')_fe('%c% var_fe %c% ')_c('%c% var_c %c%')_s('%c% var_s %c%')_sm(' %c% var_spm %c% ')' %c% '_trl(' %c% var_tr %c% ')_l(' %c% specl %c% ')_' %c% scriptv %c% spec_suf
            print(spec_ifes)
            mod <- 'DDiD_' %c% spec_ifes %c% '.Rqs'
            models <- list.files(path = cdd['p_rm'] %c% '/Clarity' %c% cfolder,pattern='\\.Rqs')
            filem <- cdd['p_rm_cl']  %c% cfolder %c% '/' %c% mod
            folder_rt <- cdd['p_rt'] %c% '/Clarity/' %c% cfolder %c% spec_ifes
            folder_rf <- cdd['p_rf'] %c% '/Clarity/' %c% cfolder %c% spec_ifes
            folder_rta <- cdd['p_rt'] %c% '/Clarity' %c% cfolder %c% '/All'
            folder_rfa <- cdd['p_rf'] %c% '/Clarity' %c% cfolder %c% '/All'
            
            var_y <- unlist(spec_y)
            var_y_missing <- var_y[!(var_y %in% names(ctfpt))]
            if (length(var_y_missing)>0){ print(var_y_missing) }
            var_y <- var_y[var_y %in% names(ctfpt)]
            var_y <- 'c(' %c% paste0(var_y,collapse=', ') %c% ') '
            # print(3)
            if ('est' %in% spec_task){
              res <- my_feols(
                                df=ctfpt,
                                labels=labels_sub,
                                form=list('y'=var_y,
                                          'i'=spec_i[[var_i]],
                                          'fe'=spec_fe[[var_fe]],
                                          'cluster'=spec_cluster[[var_c]],
                                          's'=spec_split[[var_s]],
                                          'spm'=var_spm),
                                filem=filem,
                                lean=spec_lean
                              )
              if (res$rc=='error'){
                print('estimation error')
                next
              }
              est <- res$est
            }
            
            if ('tbl_sum' %in% spec_task){
              tbl_all <- my_etable_df(est)
              openxlsx::write.xlsx(as_tibble(tbl_all), folder_rta %c% '/DDiD_' %c% spec_ifes %c% '.xlsx', row.names=FALSE, overwrite=TRUE)
            }
          }
        }
      }
    }
  }
}

rm(ctfpt)
rm(ctfp)
gc()
mem_used()




# ===============================================================================================================================
# ===============================================================================================================================
# SECTION 2: TRANSACTIONS ANALYSIS - EVENT APPROACH  --- Collapsed pre and post period with conditioning
# ===============================================================================================================================
# ===============================================================================================================================

###########################Begin Transaction dataset build###########################
ctfe <- read_parquet(cdd['p_d_cl_pa'] %c% cfolder %c% '/ctfe.parquet') %>% filter((typea=='C1')|(typea=='C2')|(typea=='C3')|(typea=='01'))

# Filter the data
ctfe <- ctfe %>% group_by(cohort) %>% mutate(anytreated = max(treated)) %>% ungroup
ctfe <- ctfe %>% filter(anytreated == 1)
# #Restrict observations based on the available data
ctfe <- ctfe %>% filter((open >= as.Date('2014-07-01'))&(open < as.Date('2020-01-01')))
ctfe <- ctfe %>% filter((distance<5000)&(ic>5000))


winsor <- 0.03
# Create some factor variables.
ctfe <- ctfe %>% mutate(
                        bal_highprin = case_when(bal_highprin>0 ~ bal_highprin, (bal_highprin==0) ~ NA_real_),
                        chargeoff = case_when(chargeoff>0 ~ chargeoff, (chargeoff==0) ~ NA_real_),
                        chargeoff_any = 1*(!is.na(chargeoff)),
                        paidlate = 1*(duration_e2a > 0.2),
                        
                        date=dateq,
                        typeas = case_when(typea=='C1' ~ 'Payday Int',typea=='C2' ~ 'Payday Phys',typea=='C5' ~ 'Rent 2 Own',
                                           typea=='C3' ~ 'Installment Int',typea=='C6' ~ 'Installment Phys',typea=='01' ~ 'Unsecured'),
                        typea_age = case_when(age_2g==0~typeas %c% ' x Age lte35',age_2g==1~typeas %c% ' x Age gt35'),
                        typea_g = case_when(typea=='C1' ~ 'Payday',typea=='C2' ~ 'Payday',typea=='C5' ~ 'Rent 2 Own',
                                            typea=='C3' ~ 'Installment',typea=='C6' ~ 'Installment',typea=='01' ~ 'Unsecured'),
                        d50_ic50_19=1*((distance<5000)&(ic>5000)&(dateq<as.Date('2020-01-01'))&(open<as.Date('2020-01-01'))),
                        typeag_age = case_when(age_2g==0~typea_g %c% ' x Age lte35',age_2g==1~typea_g %c% ' x Age gt35'),
                        etimey = floor(etime /4),
                        post = 1*(etime>=0),
                        state=state_y,
                        na=1
                      )
# need to winsorize
ctfe <- ctfe %>% group_by(typea) %>% mutate(
                                            #Winsorize outcome variables.
                                            bal_highprin_w = winsorize_x(bal_highprin, cut_l=0.0, cut_u=winsor),
                                            ) %>% ungroup()

ctfe <- ctfe %>% group_by(cohort,id,typea) %>% mutate(
                                                      trans_na = n(),
                                                      trans_prepost = 1*(n_distinct(post)==2),
                                                      trans_pre = 1*(min(post)==0)
                                                      ) %>% ungroup()

# Aggregate the dataset into a simple static pre-post dataset
ctfes <- ctfe %>% group_by(cohort,id,typea,typeas,post) %>% summarise(
  
                                                              bal_highprin_w_mu = mean(bal_highprin_w),
                                                              paidlate_mu = mean(paidlate),
                                                              chargeoff_mx = max(chargeoff_any),
                                                              
                                                              age_2g = first(age_2g),
                                                              pop2sqkm_wq10 = first(pop2sqkm_wq10),
                                                              state = first(state),
                                                              open = first(open),
                                                              typea_g = first(typea_g),
                                                              treated = first(treated),
                                                              typeag_age = first(typeag_age),
                                                              trans_prepost = first(trans_prepost),
                                                              trans_pre = first(trans_pre),
                                                              ) %>% ungroup()
ctfes <- ctfes %>% mutate(
                          age_2gr = 1 - age_2g,
                          payday = 1*(typea_g=='Payday'),
                          payday_pp = 1*((typea_g=='Payday')&(trans_prepost==1)),
                          payday_p = 1*((typea_g=='Payday')&(trans_pre==1)),
                        )

# Apply variable labels
labels_sub <- labels_all[names(labels_all) %in% names(ctfes)]
var_label(ctfes) <- labels_sub


###################set up regressions we want to loop through########################## 
spec_task = c('est','tbl_sum')
spec_y = list(
              #Table 5
              'terms'=c('bal_highprin_w_mu'),
              'payback2'=c('paidlate_mu','chargeoff_mx')
            )  
spec_i = list(
                'Tp' = 'treated:post'
              )
spec_fe = list(
                'ci_cppd10'='cohort^id + cohort^post^pop2sqkm_wq10'
                )

#Dictionary to handle the splits by defining the variable and value labels and the base category (for marginal effects)
spec_split = list(
                  'age_2g_b1'=list('var'='age_2g','var_label'='Age','base'=1,'val_label'=c('Age lte35'=0,'Age gt35'=1)),
                  'na'='na'
                ) 
spec_split_m = c('na','s','ig')  #'na','s','is','ig'
interaction_types = list('is'='specific','ig'='generic')
spec_treat_rest <- c(
                      'payday',
                      'payday_pp',
                      'payday_p'
                    )
spec_cluster <- list(
                      'coh'='cohort'
                    )
spec_lean <- TRUE  #TRUE means small files or FALSE for large files. 
spec_suf <- '_ds(ctfess)'
gc()
mem_used()



##############run regressions specified above and save output################
for (var_tr in spec_treat_rest) {
  ctfest <- ctfes %>% filter(get(var_tr)==1)
  for (var_s in names(spec_split)){
    for (var_spm in spec_split_m) {
      for (var_i in names(spec_i)){
        for (var_c in names(spec_cluster)){
          for (var_fe in names(spec_fe)){
            if ((var_s == 'na')&(var_spm!='na')){ next }
            if ((var_s != 'na')&(var_spm=='na')){ next }
            # Create folder for model, tables, figures
            specl <- ifelse(spec_lean, 'L', 'F')
            spec_ifes <- 'i(' %c% var_i %c% ')_fe('%c% var_fe %c% ')_c('%c% var_c %c%')_s('%c% var_s %c%')_sm(' %c% var_spm %c% ')' %c% '_trl(' %c% var_tr %c% ')_l(' %c% specl %c% ')_' %c% scriptv %c% spec_suf
            print(spec_ifes)
            mod <- 'DDiD_' %c% spec_ifes %c% '.Rqs'
            models <- list.files(path = cdd['p_rm'] %c% '/Clarity' %c% cfolder,pattern='\\.Rqs')
            filem <- cdd['p_rm_cl']  %c% cfolder %c% '/' %c% mod
            folder_rt <- cdd['p_rt'] %c% '/Clarity/' %c% cfolder %c% spec_ifes
            folder_rf <- cdd['p_rf'] %c% '/Clarity/' %c% cfolder %c% spec_ifes
            folder_rta <- cdd['p_rt'] %c% '/Clarity' %c% cfolder %c% '/All'
            folder_rfa <- cdd['p_rf'] %c% '/Clarity' %c% cfolder %c% '/All'
            
            var_y <- unlist(spec_y)
            var_y_missing <- var_y[!(var_y %in% names(ctfest))]
            if (length(var_y_missing)>0){ print(var_y_missing) }
            var_y <- var_y[var_y %in% names(ctfest)]
            var_y <- 'c(' %c% paste0(var_y,collapse=', ') %c% ') '
            # print(3)
            if ('est' %in% spec_task){
              res <- my_feols(
                df=ctfest,
                labels=labels_sub,
                form=list('y'=var_y,
                          'i'=spec_i[[var_i]],
                          'fe'=spec_fe[[var_fe]],
                          'cluster'=spec_cluster[[var_c]],
                          's'=spec_split[[var_s]],
                          'spm'=var_spm),
                filem=filem,
                lean=spec_lean
              )
              if (res$rc=='error'){
                print('estimation error')
                next
              }
              est <- res$est
            }
            
            if ('tbl_sum' %in% spec_task){
              tbl_all <- my_etable_df(est)
              openxlsx::write.xlsx(as_tibble(tbl_all), folder_rta %c% '/DDiD_' %c% spec_ifes %c% '.xlsx', row.names=FALSE, overwrite=TRUE)
            }
          }
        }
      }
    }
  }
}


rm(ctfe)
rm(ctfes)
rm(ctfest)
gc()
mem_used()




# # ===============================================================================================================================
# # ===============================================================================================================================
# # SECTION 4: Custom Tables.
# # ===============================================================================================================================
# # ===============================================================================================================================
ofolder_t <- cdd['p_rt_cl'] %c% cfolder %c% '/Custom'

labels_all[['interq_0']] = 'Age lte35'
labels_all[['interq_1']] = 'Income gt2.5k'
labels_all[['bal_highprin_ln']] = 'Balance ln'
labels_all[['paidearly']] = 'Paid Early'
labels_all[['paidlate']] = 'Paid Late'
labels_all[['inq_npre_any']] = 'Inq Before Open'
labels_all[['temp_inter']] = 'Sample'

#Set the labels used for all tables
my.style.df = style.df(fixef.title = "", fixef.suffix = " FE", stats.title = "", slopes.title = "")
my.style.tex = style.tex(main="base", fixef.suffix = " FE", slopes.title = "",
                         tablefoot = FALSE,line.bottom="\\bottomrule",line.top="\\toprule",
                         depvar.title = "",var.title="\\midrule", fixef.title = '\\midrule', stats.title = '\\midrule')
setFixest_etable(digits = "r4", digits.stats = "r4", style.df = my.style.df, style.tex = my.style.tex, fitstat = c('n','r2'))
setFixest_dict(unlist(labels_all))


#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Table 3: Income heterogeneity in the effect of plasma openings on non-bank credit inquiries
estb <- qread(cdd['p_rm_cl'] %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(income_3g_b2_3k)_sm(s)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
my_etable(list(estb[lhs='inq_payday_any'],estb[lhs='inq_inst_any']), file = ofolder_t %c% '/Table 3.tex', replace = TRUE)
rm(estb)
gc()
mem_used()

#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Table 4: Effect of plasma openings on non-bank credit inquiries, storefront vs. Internet
# Table for other forms of interaction heterogeneity in clarity treatment effect (full sample)
esta <- qread(cdd['p_rm_cl'] %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(na)_sm(na)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
estb <- qread(cdd['p_rm_cl'] %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(age_2g_b1)_sm(ig)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
my_etable(list(esta[lhs='inq_payday_s_any'], esta[lhs='inq_payday_i_any'],
               estb[lhs='inq_payday_s_any'], estb[lhs='inq_payday_i_any'],
               esta[lhs='inq_inst_s_any'], esta[lhs='inq_inst_i_any'],
               estb[lhs='inq_inst_s_any'], estb[lhs='inq_inst_i_any']), file = ofolder_t %c% '/Table 4.tex', replace = TRUE)
rm(esta)
rm(estb)
gc()
mem_used()

#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Table 5: Effect of plasma opening on payday loan balances and delinquencies, by age
estta1 <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tp)_fe(ci_cppd10)_c(coh)_s(na)_sm(na)_trl(payday_pp)_l(L)_v8_ds(ctfess).Rqs')
estta2 <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tp)_fe(ci_cppd10)_c(coh)_s(age_2g_b1)_sm(ig)_trl(payday_pp)_l(L)_v8_ds(ctfess).Rqs')
my_etable(list(estta1[lhs='bal_highprin_w_mu'],estta2[lhs='bal_highprin_w_mu'],
               estta1[lhs='paidlate_mu'],estta2[lhs='paidlate_mu'],
               estta1[lhs='chargeoff_mx'],estta2[lhs='chargeoff_mx']), file= ofolder_t %c% '/Table 5.tex')
rm(estta1)
rm(estta2)
gc()
mem_used()

#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Table A7: Effect of plasma openings on non-bank credit inquiries, by age
esta <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(na)_sm(na)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
estb <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(age_2g_b1)_sm(ig)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
my_etable(list(esta[lhs='inq_payday_any'],esta[lhs='inq_inst_any'],
               estb[lhs='inq_payday_any'],estb[lhs='inq_inst_any']), file=ofolder_t %c% '/Table A7.tex')
rm(esta)
rm(estb)
gc()
mem_used()

#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Table A8: Effect of plasma openings on payday loan transactions, by age
# CTFP Results by age (with both split and incremental tests.
esttas <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(age_2g_b1)_sm(s)_trl(d50_ic50_19)_l(L)_v8_ds(ctfp).Rqs')
estta <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(age_2g_b1)_sm(ig)_trl(d50_ic50_19)_l(L)_v8_ds(ctfp).Rqs')
my_etable(list(esttas[sample='Age lte35',lhs='ntrans_da_payday_any'],
               esttas[sample='Age gt35',lhs='ntrans_da_payday_any'],
               estta[lhs='ntrans_da_payday_any']), file= ofolder_t %c% '/Table A8.tex')
rm(esttas)
rm(estta)
gc()
mem_used()

#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Table A9: Robustness of the effect of plasma openings on non-bank credit inquiries to alternative control groups
estb <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(na)_sm(na)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
estp <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(na)_sm(na)_trl(ct_past)_l(L)_v8_ds(cifp).Rqs')
estf <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(na)_sm(na)_trl(ct_future)_l(L)_v8_ds(cifp).Rqs')
estpf <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(na)_sm(na)_trl(ct_fp2)_l(L)_v8_ds(cifp).Rqs')
my_etable(list(estb[lhs='inq_payday_any'],estb[lhs='inq_inst_any'],
               estp[lhs='inq_payday_any'],estp[lhs='inq_inst_any'],
               estf[lhs='inq_payday_any'],estf[lhs='inq_inst_any'],
               estpf[lhs='inq_payday_any'],estpf[lhs='inq_inst_any']), file=ofolder_t %c% '/Table A9.tex')
rm(estb)
rm(estp)
rm(estf)
rm(estpf)
gc()
mem_used()

#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Table A10: Other forms of heterogeneity in the effect of plasma openings on non-bank credit inquiries
esta <- qread(cdd['p_rm_cl'] %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(num_inq_npre_gt1)_sm(ig)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
estb <- qread(cdd['p_rm_cl'] %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(banked_2g_low)_sm(ig)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
estc <- qread(cdd['p_rm_cl'] %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(unemprt_4g_top75)_sm(ig)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
my_etable(list(esta[lhs='inq_payday_any'],
               estb[lhs='inq_payday_any'],
               estc[lhs='inq_payday_any'],
               esta[lhs='inq_inst_any'],
               estb[lhs='inq_inst_any'],
               estc[lhs='inq_inst_any']), file = ofolder_t %c% '/Table A10.tex', replace = TRUE)
rm(esta)
rm(estb)
rm(estc)
gc()
mem_used()

#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Table A11: Density heterogeneity in the effect of plasma openings on non-bank credit inquiries
estz <- qread(cdd['p_rm_cl'] %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(density_3g_low)_sm(s)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
my_etable(list(estz[lhs='inq_payday_any'],estz[lhs='inq_inst_any']), file = ofolder_t %c% '/Table A11.tex', replace = TRUE)
rm(estz)
gc()
mem_used()

#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Appendix C test of impact of selection on the dependent variable:
# Table C1: Test of sample selection on the effect of plasma openings on non-bank credit inquiries
estz <- qread(cdd['p_rm_cl'] %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(inq_npre_any)_sm(ig)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
my_etable(list(estz[lhs='inq_payday_any'],estz[lhs='inq_inst_any']), file = ofolder_t %c% '/Table C1.tex', replace = TRUE)
rm(estz)
gc()
mem_used()


# # ===============================================================================================================================
# # ===============================================================================================================================
# # SECTION 5: Graphs.
# # ===============================================================================================================================
# # ===============================================================================================================================
ofolder_f <- cdd['p_rf_cl'] %c% cfolder %c% '/Custom'


labels_all[['interq_1']] = 'Income gt2.5k'
labels_all[['bal_highprin_ln']] = 'Balance ln'
labels_all[['paidearly']] = 'Paid Early'
labels_all[['paidlate']] = 'Paid Late'
labels_all[['inq_npre_any']] = 'Inq Before Open'
labels_all[['temp_inter']] = 'Sample'

mscale = c('inq_payday_any'=0.034,'inq_inst_any'=0.047,'inq_payday_gt1'=0.022/0.034,'inq_inst_gt1'=0.028/0.047)
mscale_ctfp = c('ntrans_da_payday_any' = 0.0476)


#Set the labels used for all tables
my.style.df = style.df(fixef.title = "", fixef.suffix = " FE", stats.title = "", slopes.title = "")
my.style.tex = style.tex(main="base", fixef.suffix = " FE", slopes.title = "",
                         tablefoot = FALSE,line.bottom="\\bottomrule",line.top="\\toprule",
                         depvar.title = "",var.title="\\midrule", fixef.title = '\\midrule', stats.title = '\\midrule')
setFixest_etable(digits = "r4", digits.stats = "r4", style.df = my.style.df, style.tex = my.style.tex, fitstat = c('n','r2'))
setFixest_dict(unlist(labels_all))


#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Figure 4: Effect of plasma openings on non-bank credit inquiries
# Figure 4 - Panel A:
#CIFP__Baseline Results
est <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(na)_sm(na)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
# est_ext <- est[lhs=c('inq_payday_any','inq_inst_any','inq_payday_gt1$','inq_inst_gt1$')]
object <- est[lhs=c('inq_payday_any','inq_inst_any')]
gg <- my_iplot(est[lhs=c('inq_payday_any','inq_inst_any')],labels=labels_all,rescale=mscale,
               file=ofolder_f %c% '/Figure 4: Panel A: DDiDs_cifp_base_gt0.png')
#CIFP__Het by Age Results
labels_all[['interq_0']] = 'Treated x Age lte35'
esta <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/Figure 4: Panel A: DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(age_2g_b1)_sm(s)_trl(d50_ic50_19)_l(L)_v8_ds(cifp).Rqs')
# Figure 4 - Panel B:
gg <- my_iplot(esta[lhs=c('inq_payday_any')],labels=labels_all,split='s',rescale=mscale,
               file=ofolder_f %c% '/Figure 4: Panel B: DDiDs_sp_cifp_payday_age.png')
# Figure 4 - Panel C:
gg <- my_iplot(esta[lhs=c('inq_inst_any')],labels=labels_all,split='s',rescale=mscale,
               file=ofolder_f %c% '/Figure 4: Panel C: DDiDs_sp_cifp_inst_age.png')


#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Figure 5: Effect of plasma opening on payday transactions, by age
# CTFP Results by Age
labels_all[['interq_0']] = 'Treated x Age lte35'
estta <- qread(cdd['p_rm_cl'] %c% cfolder %c% '/DDiD_i(Tty)_fe(ci_ctpd10_sd)_c(coh)_s(age_2g_b1)_sm(s)_trl(d50_ic50_19)_l(L)_v8_ds(ctfp).Rqs')
gg <- my_iplot(estta[lhs=c('ntrans_da_payday_any')],labels=labels_all,split='s',rescale=mscale_ctfp,
               file=ofolder_f %c% '/Figure 5: DDiDs_sp_ctfp_payday_age.png')




